########
# Predicted probabilities of different agreement profiles
#######

#
rm(list=ls())


# load libraries
packs = c('lmtest', 'tidyverse', 'hrbrthemes', 'stringr', 'foreign', 
          'ggrepel', "here")

#
source("r/LoadPkg.R")
#
loadPkg(packs)

# read data
load('Data/conjoint-nodup.rda')
df = conjoint_nodup
rm(conjoint_nodup)




# High-authoritarian pred probs -------------------------------------------

# model formula 
covars_fac = c('tierras_fac*rearing_agg', 'elections_fac*rearing_agg',
               'reint_fac*rearing_agg', 'just_fac*rearing_agg', 
               'repar_fac*rearing_agg', 'drugs_fac*rearing_agg')

form_fac = formula(paste('chosen_dummy ~ ', 
                         paste(covars_fac, collapse = '+')))


# get pred probs
mod_fac = 
  lm(form_fac, data = df)
df$fitted = predict(mod_fac, type = 'response')


# scenario: actual agreement
newdata_agreement = data.frame(tierras_fac = 'land-medium', 
                       elections_fac = 'elections-high', 
                       reint_fac = 'reinter-medium', 
                       just_fac = 'justice-low', 
                       repar_fac ='repar-medium', 
                       drugs_fac = 'drug-air', 
                       rearing_agg = 4)



# pred probs for actual agreement
actual_preds = predict(mod_fac, 
                       newdata = newdata_agreement, 
                       interval = 'confidence')



# variations on actual agreement
## better land
land = newdata_agreement
land$tierras_fac = 'land-high'
impL = predict(mod_fac, newdata = land, interval = 'confidence')

## less impunity
imp = newdata_agreement
imp$just_fac = 'justice-high'
impJ = predict(mod_fac, newdata = imp, interval = 'confidence')

## less electoral
elec = newdata_agreement
elec$elections_fac = 'elections-low'
impElec = predict(mod_fac, newdata = elec, interval = 'confidence')

## less reintegration
reinter = newdata_agreement
reinter$reint_fac = 'reinter-low'
impRe = predict(mod_fac, newdata = reinter, interval = 'confidence')

## reparations high
repar = newdata_agreement
repar$repar_fac = 'repar-high'
impRest = predict(mod_fac, newdata = repar, interval = 'confidence')

## softer drug policy
drug = newdata_agreement
drug$drugs_fac = 'drug-sub'
impDr = predict(mod_fac, newdata = drug, interval = 'confidence')


## plot data
impDat = rbind(actual_preds, impL, impJ, impElec, 
               impDr, impRest, impRe) %>% 
  as.data.frame()
impDat$var = c('Actual agreement',
               'Larger land transfer', 
               'All rebels face jail time',
               'No elections', 
               'Drug substitution', 
               'Monetary reparations', 
               'No economic aid')

# plotting colors
impDat$colors = ifelse(impDat$var == 'Actual agreement', 'yes', 'no')

# separate object for facet plot
impDat2 = impDat




# Low authoritarian pred probs --------------------------------------------

# actual agreement
newdata_agreement = data.frame(tierras_fac = 'land-medium', 
                               elections_fac = 'elections-high', 
                               reint_fac = 'reinter-medium', 
                               just_fac = 'justice-low', 
                               repar_fac ='repar-medium', 
                               drugs_fac = 'drug-air', 
                               rearing_agg = 1)



# actual agreement
actual_preds = predict(mod_fac, 
                       newdata = newdata_agreement, 
                       interval = 'confidence')


# variations on actual agreement
## better land
land = newdata_agreement
land$tierras_fac = 'land-high'
impL = predict(mod_fac, newdata = land, interval = 'confidence')

## impunity
imp = newdata_agreement
imp$just_fac = 'justice-high'
impJ = predict(mod_fac, newdata = imp, interval = 'confidence')

## electoral
elec = newdata_agreement
elec$elections_fac = 'elections-low'
impElec = predict(mod_fac, newdata = elec, interval = 'confidence')

## reintegration
reinter = newdata_agreement
reinter$reint_fac = 'reinter-low'
impRe = predict(mod_fac, newdata = reinter, interval = 'confidence')

## reparations high
repar = newdata_agreement
repar$repar_fac = 'repar-high'
impRest = predict(mod_fac, newdata = repar, interval = 'confidence')

drug = newdata_agreement
drug$drugs_fac = 'drug-sub'
impDr = predict(mod_fac, newdata = drug, interval = 'confidence')


## plot data
impDat = rbind(actual_preds, impL, impJ, impElec, 
               impDr, impRest, impRe) %>% 
  as.data.frame()
impDat$var = c('Actual agreement',
               'Larger land transfer', 
               'All rebels face jail time',
               'No elections', 
               'Drug substitution', 
               'Monetary reparations', 
               'No economic aid')

# plotting colors
impDat$colors = ifelse(impDat$var == 'Actual agreement', 'yes', 'no')


# facet plot
impDat$facet = 'Low child-rearing score'
impDat2$facet = 'High child-rearing score'

pDat = rbind(impDat, impDat2)

#
ggplot(pDat) + geom_hline(yintercept = .5,
           colour = gray(1/2),
           lty = 2) +
  geom_pointrange(aes(x = var, y = fit, ymin = lwr, ymax = upr,
                      fill = colors, color = colors, shape = colors),
                  lwd = 1/2, fatten = 6) +
  coord_flip() +
  theme_ipsum_rc(grid = 'Y') +
  scale_y_percent() +
  scale_fill_manual(values = c('#00A08A', '#FF0000')) +
  scale_color_manual(values = c('#00A08A', '#FF0000')) +
  theme(legend.position = 'none', strip.text = element_text(face = 'bold')) +
  labs(y = 'Predicted probability of agreement being chosen',
       x = '') +
  facet_wrap(vars(facet), nrow = 2)

ggsave(here("repFile", "paper", "figures", "pred-prob-facet.pdf"), 
       device = cairo_pdf)
ggsave(here("repFile", "JOP submission", "pred-prob-facet.eps"), 
       device = cairo_ps)
ggsave(here("repFile", "JOP submission", "pred-prob-facet.pdf"), 
       device = cairo_pdf)


# Comparing voter types --------------------------------------------------------

covars_fac = c('tierras_fac', 'elections_fac', 'reint_fac', 'just_fac', 
               'repar_fac', 'drugs_fac')

form_fac = formula(paste('chosen_dummy ~ ', paste(covars_fac, collapse = '+')))


# Predicted Probabilities for subsets fo voters
mod_yes = 
  lm(form_fac, data = df[df$vote_choice == 'A favor',])

mod_no = 
  lm(form_fac, data = df[df$vote_choice == 'En contra',])

mod_abs = 
  lm(form_fac, data = df[df$vote_choice == 'Did not vote',])



## improvements function
improve_agreements = function(model) 
{
  
  # actual agreement
  newdata_agreement = data.frame(tierras_fac = 'land-medium', 
                                 elections_fac = 'elections-high', 
                                 reint_fac = 'reinter-medium', 
                                 just_fac = 'justice-low', 
                                 repar_fac ='repar-medium', 
                                 drugs_fac = 'drug-air')
  
  # actual agreement
  actual_preds = predict(model, 
                         newdata = newdata_agreement, 
                         interval = 'confidence')
  
  ## better land
  land = newdata_agreement
  land$tierras_fac = 'land-high'
  impL = predict(model, newdata = land, interval = 'confidence')
  
  ## impunity
  imp = newdata_agreement
  imp$just_fac = 'justice-high'
  impJ = predict(model, newdata = imp, interval = 'confidence')
  
  ## electoral
  elec = newdata_agreement
  elec$elections_fac = 'elections-low'
  impElec = predict(model, newdata = elec, interval = 'confidence')
  
  ## reintegration
  reinter = newdata_agreement
  reinter$reint_fac = 'reinter-low'
  impRe = predict(model, newdata = reinter, interval = 'confidence')
  
  ## reparations high
  repar = newdata_agreement
  repar$repar_fac = 'repar-high'
  impRest = predict(model, newdata = repar, interval = 'confidence')
  
  drug = newdata_agreement
  drug$drugs_fac = 'drug-sub'
  impDr = predict(model, newdata = drug, interval = 'confidence')
  
  
  ## plot data
  impDat = rbind(actual_preds, impL, impJ, impElec, 
                 impDr, impRest, impRe) %>% 
    as.data.frame()
  impDat$var = c('Actual agreement',
                 'Larger land transfer', 
                 'All rebels face jail time',
                 'No elections', 
                 'Drug substitution', 
                 'Monetary reparations', 
                 'No economic aid')
  
  # plotting colors
  impDat$colors = ifelse(impDat$var == 'Actual agreement', 'yes', 'no')
  
  
  return(impDat)
  
}

# get improved agreements
imp_yes = 
  improve_agreements(mod_yes) %>% 
  mutate(type = 'Voted Yes')
imp_no = 
  improve_agreements(mod_no) %>% 
  mutate(type = 'Voted No')
imp_abs = 
  improve_agreements(mod_abs) %>% 
  mutate(type = 'Abstained')

# combine for plotting
pDat = rbind(imp_yes, imp_no, imp_abs)

# reorder
pDat$type = factor(pDat$type, levels = c('Voted Yes', 'Voted No', 
                                         'Abstained'))

# plot
pDat %>% 
  ggplot() +
  geom_hline(yintercept = .5,
             colour = gray(1/2),
             lty = 2) +
  geom_pointrange(aes(x = var, y = fit, ymin = lwr, ymax = upr,
                      fill = colors, color = colors, shape = colors),
                  lwd = 1/2, fatten = 6) +
  coord_flip() +
  theme_ipsum_rc(grid = 'Y') +
  scale_y_percent() +
  scale_fill_manual(values = c('#00A08A', '#FF0000')) +
  scale_color_manual(values = c('#00A08A', '#FF0000')) +
  theme(legend.position = 'none', strip.text = element_text(face = 'bold')) +
  labs(y = 'Predicted probability of agreement being chosen',
       x = '') +
  facet_wrap(vars(type), nrow = 3)

ggsave(here("repFile", "paper", "figures", "improvement-vote.pdf"),
       device = cairo_pdf)
ggsave(here("repFile", "JOP submission", "improvement-vote.eps"),
       device = cairo_ps)
ggsave(here("repFile", "JOP submission", "improvement-vote.pdf"),
       device = cairo_pdf)